home *** CD-ROM | disk | FTP | other *** search
/ Chip 2007 January, February, March & April / Chip-Cover-CD-2007-02.iso / Pakiet multimedia / Grafika i zdjecia / Edytory grafiki rastrowej i wektorowej / Inscape 0.44.1 / Inkscape-0.44.1-1.win32.exe / share / extensions / bezmisc.py < prev    next >
Text File  |  2006-09-06  |  8KB  |  251 lines
  1. #!/usr/bin/env python
  2. '''
  3. Copyright (C) 2005 Aaron Spike, aaron@ekips.org
  4.  
  5. This program is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU General Public License as published by
  7. the Free Software Foundation; either version 2 of the License, or
  8. (at your option) any later version.
  9.  
  10. This program is distributed in the hope that it will be useful,
  11. but WITHOUT ANY WARRANTY; without even the implied warranty of
  12. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13. GNU General Public License for more details.
  14.  
  15. You should have received a copy of the GNU General Public License
  16. along with this program; if not, write to the Free Software
  17. Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18. '''
  19.  
  20. import math, cmath
  21.  
  22. def rootWrapper(a,b,c,d):
  23.     if a:
  24.         #TODO: find a new cubic solver and put it here
  25.           #return solveCubicMonic(b/a,c/a,d/a)
  26.         return ()
  27.     elif b:
  28.         det=c**2.0-4.0*b*d
  29.         if det:
  30.             return (-c+cmath.sqrt(det))/(2.0*b),(-c-cmath.sqrt(det))/(2.0*b)
  31.         else:
  32.             return -c/(2.0*b),
  33.     elif c:
  34.         return 1.0*(-d/c),
  35.     return ()
  36.  
  37. def bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3))):
  38.     #parametric bezier
  39.     x0=bx0
  40.     y0=by0
  41.     cx=3*(bx1-x0)
  42.     bx=3*(bx2-bx1)-cx
  43.     ax=bx3-x0-cx-bx
  44.     cy=3*(by1-y0)
  45.     by=3*(by2-by1)-cy
  46.     ay=by3-y0-cy-by
  47.  
  48.     return ax,ay,bx,by,cx,cy,x0,y0
  49.     #ax,ay,bx,by,cx,cy,x0,y0=bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)))
  50.  
  51. def linebezierintersect(((lx1,ly1),(lx2,ly2)),((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3))):
  52.     #parametric line
  53.     dd=lx1
  54.     cc=lx2-lx1
  55.     bb=ly1
  56.     aa=ly2-ly1
  57.  
  58.     if aa:
  59.         coef1=cc/aa
  60.         coef2=1
  61.     else:
  62.         coef1=1
  63.         coef2=aa/cc
  64.  
  65.     ax,ay,bx,by,cx,cy,x0,y0=bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)))
  66.     #cubic intersection coefficients
  67.     a=coef1*ay-coef2*ax
  68.     b=coef1*by-coef2*bx
  69.     c=coef1*cy-coef2*cx
  70.     d=coef1*(y0-bb)-coef2*(x0-dd)
  71.  
  72.     roots = rootWrapper(a,b,c,d)
  73.     retval = []
  74.     for i in roots:
  75.         if type(i) is complex and i.imag==0:
  76.             i = i.real
  77.         if type(i) is not complex and 0<=i<=1:
  78.             retval.append(i)
  79.     return retval
  80.  
  81. def bezierpointatt(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)),t):
  82.     ax,ay,bx,by,cx,cy,x0,y0=bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)))
  83.     x=ax*(t**3)+bx*(t**2)+cx*t+x0
  84.     y=ay*(t**3)+by*(t**2)+cy*t+y0
  85.     return x,y
  86.  
  87. def bezierslopeatt(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)),t):
  88.     ax,ay,bx,by,cx,cy,x0,y0=bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)))
  89.     dx=3*ax*(t**2)+2*bx*t+cx
  90.     dy=3*ay*(t**2)+2*by*t+cy
  91.     return dx,dy
  92.  
  93. def beziertatslope(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)),(dy,dx)):
  94.     ax,ay,bx,by,cx,cy,x0,y0=bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)))
  95.     #quadratic coefficents of slope formula
  96.     if dx:
  97.         slope = 1.0*(dy/dx)
  98.         a=3*ay-3*ax*slope
  99.         b=2*by-2*bx*slope
  100.         c=cy-cx*slope
  101.     elif dy:
  102.         slope = 1.0*(dx/dy)
  103.         a=3*ax-3*ay*slope
  104.         b=2*bx-2*by*slope
  105.         c=cx-cy*slope
  106.     else:
  107.         return []
  108.  
  109.     roots = rootWrapper(0,a,b,c)
  110.     retval = []
  111.     for i in roots:
  112.         if type(i) is complex and i.imag==0:
  113.             i = i.real
  114.         if type(i) is not complex and 0<=i<=1:
  115.             retval.append(i)
  116.     return retval
  117.  
  118. def tpoint((x1,y1),(x2,y2),t):
  119.     return x1+t*(x2-x1),y1+t*(y2-y1)
  120. def beziersplitatt(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)),t):
  121.     m1=tpoint((bx0,by0),(bx1,by1),t)
  122.     m2=tpoint((bx1,by1),(bx2,by2),t)
  123.     m3=tpoint((bx2,by2),(bx3,by3),t)
  124.     m4=tpoint(m1,m2,t)
  125.     m5=tpoint(m2,m3,t)
  126.     m=tpoint(m4,m5,t)
  127.     
  128.     return ((bx0,by0),m1,m4,m),(m,m5,m3,(bx3,by3))
  129.  
  130. '''
  131. Approximating the arc length of a bezier curve
  132. according to <http://www.cit.gu.edu.au/~anthony/info/graphics/bezier.curves>
  133.  
  134. if:
  135.     L1 = |P0 P1| +|P1 P2| +|P2 P3| 
  136.     L0 = |P0 P3|
  137. then: 
  138.     L = 1/2*L0 + 1/2*L1
  139.     ERR = L1-L0
  140. ERR approaches 0 as the number of subdivisions (m) increases
  141.     2^-4m
  142.  
  143. Reference:
  144. Jens Gravesen <gravesen@mat.dth.dk>
  145. "Adaptive subdivision and the length of Bezier curves"
  146. mat-report no. 1992-10, Mathematical Institute, The Technical
  147. University of Denmark. 
  148. '''
  149. def pointdistance((x1,y1),(x2,y2)):
  150.     return math.sqrt(((x2 - x1) ** 2) + ((y2 - y1) ** 2))
  151. def Gravesen_addifclose(b, len, error = 0.001):
  152.     box = 0
  153.     for i in range(1,4):
  154.         box += pointdistance(b[i-1], b[i])
  155.     chord = pointdistance(b[0], b[3])
  156.     if (box - chord) > error:
  157.         first, second = beziersplitatt(b, 0.5)
  158.         Gravesen_addifclose(first, len, error)
  159.         Gravesen_addifclose(second, len, error)
  160.     else:
  161.         len[0] += (box / 2.0) + (chord / 2.0)
  162. def bezierlengthGravesen(b, error = 0.001):
  163.     len = [0]
  164.     Gravesen_addifclose(b, len, error)
  165.     return len[0]
  166.  
  167. # balf = Bezier Arc Length Function
  168. balfax,balfbx,balfcx,balfay,balfby,balfcy = 0,0,0,0,0,0
  169. def balf(t):
  170.     retval = (balfax*(t**2) + balfbx*t + balfcx)**2 + (balfay*(t**2) + balfby*t + balfcy)**2
  171.     return math.sqrt(retval)
  172.  
  173. def Simpson(f, a, b, n_limit, tolerance):
  174.     n = 2
  175.     multiplier = (b - a)/6.0
  176.     endsum = f(a) + f(b)
  177.     interval = (b - a)/2.0
  178.     asum = 0.0
  179.     bsum = f(a + interval)
  180.     est1 = multiplier * (endsum + (2.0 * asum) + (4.0 * bsum))
  181.     est0 = 2.0 * est1
  182.     #print multiplier, endsum, interval, asum, bsum, est1, est0
  183.     while n < n_limit and abs(est1 - est0) > tolerance:
  184.         n *= 2
  185.         multiplier /= 2.0
  186.         interval /= 2.0
  187.         asum += bsum
  188.         bsum = 0.0
  189.         est0 = est1
  190.         for i in xrange(1, n, 2):
  191.             bsum += f(a + (i * interval))
  192.             est1 = multiplier * (endsum + (2.0 * asum) + (4.0 * bsum))
  193.     #print multiplier, endsum, interval, asum, bsum, est1, est0
  194.     return est1
  195.  
  196. def bezierlengthSimpson(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)), tolerance = 0.001):
  197.     global balfax,balfbx,balfcx,balfay,balfby,balfcy
  198.     ax,ay,bx,by,cx,cy,x0,y0=bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)))
  199.     balfax,balfbx,balfcx,balfay,balfby,balfcy = 3*ax,2*bx,cx,3*ay,2*by,cy
  200.     return Simpson(balf, 0.0, 1.0, 4096, tolerance)
  201.  
  202. def beziertatlength(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)), l = 0.5, tolerance = 0.001):
  203.     global balfax,balfbx,balfcx,balfay,balfby,balfcy
  204.     ax,ay,bx,by,cx,cy,x0,y0=bezierparameterize(((bx0,by0),(bx1,by1),(bx2,by2),(bx3,by3)))
  205.     balfax,balfbx,balfcx,balfay,balfby,balfcy = 3*ax,2*bx,cx,3*ay,2*by,cy
  206.     t = 1.0
  207.     tdiv = t
  208.     curlen = Simpson(balf, 0.0, t, 4096, tolerance)
  209.     targetlen = l * curlen
  210.     diff = curlen - targetlen
  211.     while abs(diff) > tolerance:
  212.         tdiv /= 2.0
  213.         if diff < 0:
  214.             t += tdiv
  215.         else:
  216.             t -= tdiv            
  217.         curlen = Simpson(balf, 0.0, t, 4096, tolerance)
  218.         diff = curlen - targetlen
  219.     return t
  220.  
  221. #default bezier length method
  222. bezierlength = bezierlengthSimpson
  223.  
  224. if __name__ == '__main__':
  225.     import timing
  226.     #print linebezierintersect(((,),(,)),((,),(,),(,),(,)))
  227.     #print linebezierintersect(((0,1),(0,-1)),((-1,0),(-.5,0),(.5,0),(1,0)))
  228.     tol = 0.00000001
  229.     curves = [((0,0),(1,5),(4,5),(5,5)),
  230.             ((0,0),(0,0),(5,0),(10,0)),
  231.             ((0,0),(0,0),(5,1),(10,0)),
  232.             ((-10,0),(0,0),(10,0),(10,10)),
  233.             ((15,10),(0,0),(10,0),(-5,10))]
  234.     '''
  235.     for curve in curves:
  236.         timing.start()
  237.         g = bezierlengthGravesen(curve,tol)
  238.         timing.finish()
  239.         gt = timing.micro()
  240.  
  241.         timing.start()
  242.         s = bezierlengthSimpson(curve,tol)
  243.         timing.finish()
  244.         st = timing.micro()
  245.  
  246.         print g, gt
  247.         print s, st
  248.     '''
  249.     for curve in curves:
  250.         print beziertatlength(curve,0.5)
  251.